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SUMMARY 


A sample heat transfer analysis is shown which includes the heat of 
fusion. The method can be used to analyze a system with nonconstant specific 
heat. The enthalpy is Introduced as an independent degree of freedom at each 
node. The user input consists of a curve of temperature as a function of 
enthalpy, which may Include a constant temperature phase change. The basic 
NASTRAN heat transfer capability is used to model the effects of latent heat 
with existing direct matrix output (DMI) and nonlinear load (N0LIN) data 
cards. Although some user care is required, the numerical stability of the 
Integration is quite good when the given recommendations are followed. The 
theoretical equations used and the NASTRAN techniques are shown in the paper. 


INTRODUCTION 


The problem of heat transfer with latent heat or nonconstant specific 
heat is one of current Interest. Methods based upon introduction of 
enthalpy^^^ and upon a very high specific heat in a small temperature range'^^ 
have been published. These methods can be implemented in either finite 
difference or finite element formulations. The numerical stability of the 
transient integration must be analyzed, since failure in this area makes the 
method very expensive, due to excessive time steps required. 
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SYMBOLS 


Values are given in SI Units* 

O 

A area, 

B capacity matrix, J/°C 

c specific heat, J/kg°C 

h enthalpy, J/kg 

K conduction matrix, J/s°C 

L latent heat, h. - h , J/kg 

JL S 

M mass, kg 

N nonlinear function 

Ng number of finite elements 

P thermal load , J/s 

T temperature, °C 

t time, s 

X coordinate, m 

p Newraark beta parameter 

Y stability parameter 

K thermal conductivity, J/s m °C 

\ eigenvalue, growth factor 

p mass density, k/m^ 

Subscripts : 
i llquldus 

s solidus 

w wall 

n time step 
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MATHEMATICAL ANALYSIS 


The derivation will be made for one-dimensional heat conduction, however 
the method is general. The temperature distribution is found from solution of 
the diffusion equation 


d^T 5T 

K =• - pc 

ax at 


( 1 ) 


When this equation is analyzed by finite element techniques, it is 
written 

[b]{t} + [k]{t} - {P} (2) 


where the terms of the B (capacitances) matrix are pcAAx and the terms of the 
K (finite element conductivities) matrix ara Ak/Ax. 

In the general case the conductivity and the heat capacity are not 
constant. In the case of phase change at constant temperature (see fig. 1) 
the method falls since the specific heat is effectively infinite. 
Introduction of the enthalpy gives two simultaneous equations 

[M]{h} + [k]{t} = {P} (3) 

T = N (h) (4) 


A Newmark beta numerical integration scheme is used, where the velocity 
terras are replaced by 


^ - (Vl - '■nVit (5) 

and the constant terras by 

T = P T^^^ + (1-p) T^ (6) 

This pararaeter p is a stability parameter. If p * 0, the method is called 
foreward differencing, and p =• 1 corresponds to backward differencing. For p 
greater than 0.5, integration of equation (2) is numerically stable for any 
mesh size and time step. Integration of (3) and (4) shows that there is a 
tendency to instability at large time steps. 
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STABILITY ANALYSIS 


The basic method of stabilizing the equations is to evaluate the terms at 
the advanced time rather than the time (t^^) . We replace equation (4) 
by ” 


- "''min - 


- h/c 


rain 


(7) 


where is the mlnimtam specific heat> The terms on the left side of (7) 

are evaluated as shown in (6) , while the terras on the right are evaluated at 
An entirely equivalent method to (7) is 

T - (yAt/c^^^) h = N (h) (8) 


If (8) is derived from (7), the parameter y would be P; however, we shall treat 
Y as an Independent stability parameter. 

Stability analysis requires a long derivation and will not be done 
here. The two basic steps are linearization and modal analysis. The 
equations are linearized by replacing the nonlinear curve by a linear fit. 
The mode analysis replaces a multi-degree-of-freedom problem with a series of 
two - degree - of - freedom problems. The short wave length modes are most 
unstable. Introduce a growth factor 


n+1 


\ T 


n 


(9) 


h 


n+1 


\ h 

n 


( 10 ) 


There are two roots to the characteristic equation 


X “ 


1-6 


-1 


and 


1 - 


pcL*- 

4N^icAt 

e 


+ Y 


"min 


-1 


( 11 ) 


Since -1 < \ < +1 for stability, 


6 > 1/2 


( 12 ) 
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and 


+ Y > 1/2 (13) 

4N <At c . 
e min 

Thus P “ Y " 0.55 gives good stability for any At. For At < pcL^/2N^K, the 
solution Is stable without the need of y» l>ut this At Is usually much Iraaller 
than needed for accuracy. 


NASTRAN RESULTS 


The analytic solution given In reference (1) was chosen, since an exact 
solution allows analysis of accuracy. This Is a one-dimensional model. The 
Initial condition Is Ice at freezing temperature. Starting at time zero, the 
end is heated to a constant temperature. The end temperature Is chosen so the 
enthalpy change after melting has the same value as the latent heat. The 
melting proceeds until such time that the water-ice Interface has moved 1.24 
meters, and then the temperature profile is examined. The model Is 1.3 meters 
long and has either 26 or 130 elements. This arrangement was chosen to allow 
comparison of results (see Table I). 


CONCLUSIONS 


The NASTRAN results are of good accuracy. By using the stability factor, 
the accuracy has not been degraded. These results can easily be adapted to 
other geometries with the finite element method. No changes were required In 
NASTRAN to solve problems with latent heat. 
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TABLE I. COMPARISON OF METHODS 


Method 

Number of Steps 


Distance 
from hot 
end , m • 


Method 

Number of Steps 

Distance 
from hot 
end f m • 


Dimensionless Temperature at the Final Time 
<T - T 3 ) / (T„ - T 3 ) 

All results are for p = y = 0*55 


Ax » 0.05 m (26 elements) 



Exact 

NASTRAN 

800 

NASTRAN 

200 

Ref 1 
800 

Ref 1 
200 

0.2 

0.8184 

0.8184 

0.8181 

0.8207 

0.8242 

0.6 

.4695 

.4694 

.4686 

.4758 

.4871 

1.2 

.0252 

.0321 

.0230 

.0335 

.0612 


Ax = 0.01 m (130 elements) 



NASTRAN 

1600 

NASTRAN 

800 

NASTRAN 

200 

Ref 1 
20000 

Ref 1 
5000 

0.2 

0.8184 

0.8184 

0.8181 

0.8189 

0.8195 

0.6 

.4694 

.4693 

.4686 

.4707 

.4726 

1.2 

.0256 

.0263 

.0273 

.0262 

.0316 
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h (enthalpy) J/kg 


Figure 1 


Latent heat; typical curve of temperature as a function 
of enthalpy, showing change of phase. The specific heat 
c = dh/dT, becomes infinite. 



